# LM_DLM_replication.R
# Replication code for
# "A Dynamic Linear Modeling Approach to Public Policy Change" by
# Matt W. Loftis and Peter B. Mortensen
# Direct inquiries to: mattwloftis@gmail.com
# R version 3.4.0
# July 2017

rm(list = ls())
library(zoo) #zoo_1.8-0
library(dlm) #dlm_1.1-4
library(data.table) #data.table_1.10.4
library(ggplot2) #ggplot2_2.2.1
source("LM_DLM_replication_helper_functions.R")


## Load replication data file
load("LM_DLM_replication_data.RData")

## Convert to time series
dat <- ts(dat, start = min(dat$year), end = max(dat$year))

## Select variables for main model
x <- dat[, c("lag_level", "chg_gdp", "lagMIP16", 
             "lag_unem", "lag_ideo", "host_lev",
             "election")]


## Define pretty labels for graphs
labs <- c("Lag Defense Spending", "Change in GDP", 
          "Lag Public Opinion", "Lag Unemployment",
          "Lag Congressional Ideology", "Hostilities", 
          "Presidential Election Year")


## Pull out dependent variable
y <- dat[,"def_change"] #annual % change in defense spending


## Define label for DV
yvar <- "Change in US Defense Spending"


## Get time periods consistent between dv & ivs
y <- na.contiguous(zoo(x = y))
x <- na.contiguous(zoo(x = x))
y <- window(y, start = start(x), end = end(x))
x <- window(x, start = start(y), end = end(y))


## Run model
results <- run.kalman(scale(y), scale(x))


## Extract labels for years included in analysis
for.labels <- as.vector(window(dat[,"year"], start = start(x), end = end(x)))


## Set up for results for graphing
t <- graph.DLM.prep(results[[1]], for.labels) #Get data to plot from smoother


## Plot each coefficient separately
sep <- split(t[[1]], t[[1]]$param) #Split output for plot by explanatory variable
stem <- "dlm_graph" #Define graph file stem
for (i in 1:length(sep)){ #Loop over explanatory variables, plotting each one
  sep[[i]]$year <- as.numeric(sep[[i]]$year) #Convert year variable to numeric for plot
  tmp <- plot.DLM.sep(coef = sep[[i]], 
                      xaxis = t[[2]], 
                      lab = labs[[i]]) #Plot individual effect coefficient over time
  ggsave(tmp, file = paste0(stem, i, ".pdf"), width = 5, height = 3) #Save graph to pdf
}



#########################################
## DIAGNOSTICS ##########################
## USING PACKAGE DLM TOOLS ##############
#########################################

## Extract standardized residuals
res <- residuals(results[[3]], sd = FALSE)


## Plot standardized residuals
plot(res, xlab = "", ylab = "Standardized residuals")


## Box test
Box.test(res, lag = 20, type = "Ljung")
sapply(1 : 20, function(i) Box.test(res, lag = i, type = "Ljung-Box")$p.value)


## Homoscedasticity of residuals
d <- 2 #number of initial variance parameters estimated
h <- round((length(res) - d) / 3)

h.stat <- sum(res[(length(res) - (h+1)):length(res)]^2) / sum(res[(d +1):(d + h)]^2)
crit.value <- qf(.975, h, h)

if(h.stat > 1) {
  h.stat < crit.value
} else {
  (1 / h.stat) < crit.value
}


## Normality
shapiro.test(as.vector(res))
